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)£> '. ABSTRACT 

T-H 

j> ■ The solution of the non-linear diffusive shock acceleration problem, where the pres- 

sure of the non-thermal population is sufficient to modify the shock hydrodynamics, is 
widely recognized as a key to understanding particle acceleration in a variety of astro- 
physical environments. We have developed a Monte Carlo technique for self-consistently 
calculating the hydrodynamic structure of oblique, steady-state shocks, together with 
the first-order Fermi acceleration process and associated non-thermal particle distribu- 
■ tions. This is the first internally consistent treatment of modified shocks that includes 

cross-field diffusion of particles. Our method overcomes the injection problem faced by 



analytic descriptions of shock acceleration, and the lack of adequate dynamic range and 



artificial suppression of cross-field diffusion faced by plasma simulations; it currently 
provides the most broad and versatile description of collisionless shocks undergoing ef- 
ficient particle acceleration. We present solutions for plasma quantities and particle 
distributions upstream and downstream of shocks, illustrating the strong differences 
observed between non-linear and test-particle cases. It is found that there are only 
marginal differences in the injection efficiency and resultant spectra for two extreme 
scattering modes, namely large-angle scattering and pitch-angle diffusion, for a wide 
range of shock parameters, i.e., for non-perpendicular subluminal shocks with field 
obliquities less than or equal to 75° and de Hoffmann- Teller frame speeds much less 
than the speed of light. 



Subject headings: Cosmic rays: general — particle acceleration — shock waves 
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1. INTRODUCTION 

The importance of shocks as generators of highly non-thermal particle distributions in helio- 
spheric and astrophysical environments has been well-documented in the literature (e.g. Axford 
1981; V61k 1984; Blandford and Eichler 1987; Jones and Ellison 1991). While direct detections of 
high energy particles are obtained via terrestrial observations of the cosmic ray flux and spacecraft 
measurements of non-thermal ions in the solar neighbourhood and in environs of planetary bow 
shocks and interplanetary travelling shocks, the existence of abundant non-thermal particle popu- 
lations in a diversity of astrophysical locales can be inferred from the prominence of non-thermal 
radiation emitted by many cosmic objects. Understanding the details of shock acceleration is of 
critical importance since many such objects emit predominantly non-thermal radiation, and indeed 
some sources are observed only because they produce non-thermal particles (e.g. radio emission 
from supernova remnants and extra-galactic jets). The first-order Fermi mechanism of diffusive 
shock acceleration is the most popular candidate for particle energization at astrophysical shocks. 
The test-particle (i.e. linear) theory (Krymsky 1977; Axford et al. 1977; Bell 1978; Blandford and 
Ostriker 1978) of this process is straightforward and yields the most important result, namely that 
a power-law with a spectral index that is relatively insensitive to the ambient conditions, is the 
natural product of collisionless shock acceleration. 

The equally important question of the efficiency of the process can only be adequately ad- 
dressed with a fully non-linear (and therefore complex) calculation. The inherent efficiency of 
shock acceleration, which is evident in observations at the Earth's bow shock (e.g. Ellison, Mobius, 
and Paschmann 1990), and in modelling of plane-parallel (e.g. see Ellison and Eichler 1985; Elli- 
son, Jones and Reynolds 1990) shocks, where the field is normal to the shock, and oblique shocks 
(e.g. Baring, Ellison and Jones 1993; Ellison, Baring and Jones 1995), implies that hydrodynamic 
feedback effects between the accelerated particles and the shock structure are very important and 
therefore essential to any complete description of the process. This has turned out to be a formidable 
task because of the wide range of spatial and energy scales that must be self-consistently included 
in a complete calculation. On the one hand, the microphysical plasma processes of the shock dis- 
sipation control injection from the thermal population and on the other hand, the highest energy 
particles (extending to at least 10 14 eV in the case of galactic cosmic rays) with extremely long 
diffusion lengths, are dynamically significant in strong shocks and feedback on the shock structure. 
Ranges of interacting scales of many orders of magnitude must be described self-consistently. 

Additional complications stem from the fact that the geometry of shocks, i.e., whether they 
are oblique or parallel, strongly affects the acceleration efficiency (e.g., Ellison, Baring, and Jones 
1995), even though the test-particle result is independent of the geometry. Observations indicate 
that interplanetary shocks, bow shocks (both planetary and from jets), the solar wind termination 
shock, and supernova remnant blast waves have a wide range of obliquities thereby rendering con- 
siderations of shock geometry salient. It turns out that the angle between the upstream magnetic 
field and the shock normal, Oeni , is a decisive parameter in determining all aspects of the shock, 
including the ability to inject and accelerate particles, and therefore has obvious observational con- 
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sequences. For instance, diffuse ions generated at quasi-parallel portions of the Earth's bow shock 
differ radically in energy content, distribution function, etc., from field aligned beams generated at 
quasi-perpendicular portions of the shock (e.g. Ipavich et al. 1988). Furthermore, the observed 
variation of radio intensity around the rim of shell-like supernova remnants may be the result of 
varying shock obliquity (e.g. Fulbright and Reynolds 1990), and the acceleration of the anomalous 
cosmic ray component at the solar wind termination shock may depend on rapid acceleration rates 
obtained in highly oblique portions of the shock (Jokipii 1992). Unfortunately, in models which 
ignore the plasma microstructure as we do here, oblique shocks are more complicated and require 
additional parameters for a complete description than do parallel (i.e. ©Bni = 0°) ones, primarily 
the degree of diffusion perpendicular to the mean ambient magnetic field direction. 

In this paper we present our method for calculating the structure of steady-state, collisionless 
shocks of arbitrary obliquity and with efficient particle injection and acceleration. The method, a 
computer simulation using Monte Carlo techniques, is an extension of our previous work on mod- 
ified parallel shocks (e.g. see Jones and Ellison 1991 and references therein), where we explored 
the properties of the non-linear modified shock scenario, and test-particle oblique shocks (Ellison, 
Baring, and Jones 1995; Baring, Ellison, and Jones 1993, 1995), where we determined the depen- 
dence of acceleration efficiency on obliquity G-Bni • These studies have been successfully applied 
to AMPTE observations near the parallel portion of the Earth's bow shock (Ellison, Mdbius, and 
Paschmann 1990), a high Mach number shock with strong modification by the accelerated ions, 
and measurements by Ulysses at highly oblique travelling interplanetary shocks in the heliosphere 
(Baring et al. 1995), which generally have low Mach numbers and therefore are well-modelled by 
linear test-particle simulations. The impressive fits obtained to the spectral data (i.e. ion distribu- 
tion functions) from each of these experiments underlines the importance of the Fermi mechanism 
and the value of the Monte Carlo technique. The present work represents the first self-consistent 
treatment of modified shocks that includes three-dimensional diffusion. 

With the Monte Carlo simulation, we self-consistently determine the average flow speed and 
magnetic field structure across the shock under the influence of accelerated particles maintaining 
constant particle, momentum, and energy fluxes at all positions from far upstream to far down- 
stream of the shock. Particles are injected upstream of the shock, propagated and diffused in the 
shock environs until they eventually leave the system. We calculate their orbits exactly as in the 
works of Decker (1988), Begelman and Kirk (1990), Ostrowski (1991) and our recent test-particle 
treatment (Ellison, Baring and Jones 1995), and make no assumption relating to the particle 
magnetic moment. Our method does not self-consistently calculate the complex plasma processes 
responsible for dissipation, but instead postulates that these processes can be adequately described 
with a simple elastic scattering relation that is assumed to be valid for all particle energies; ther- 
mal and non-thermal particles are treated identically. This simplification sacrifices the details of 
wave-particle interactions, but permits simultaneous description of the thermal plasma and the par- 
ticle injection and acceleration to the high energies associated with space plasma shocks, thereby 
satisfying the aforementioned goal of broad dynamic range. Cross-field diffusion is included via a 
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parametric description but is fully three-dimensional in contrast to hybrid or full plasma simula- 
tions with one or two ignorable dimensions which suffer from artificial suppression of cross-field 
diffusion (e.g. Jokipii, Giacalone, and Kota 1993). Simulation output includes the ion distribution 
function at all relevant positions in the shocked flow, for a range of obliquities and Mach numbers. 

Results are compared (in Section 4) for two extreme scattering modes, namely large-angle 
scattering (LAS), where the direction of a particle is isotropized in a single scattering event, and 
pitch-angle diffusion (PAD), where small changes in the angle a particle's momentum makes with 
the local magnetic field occur at each time step. The former of these extremes mimics particle 
motion in highly turbulent fields, while the latter is usually implemented in analytic treatments of 
Fermi acceleration (e.g. Kirk and Schneider 1987, but see also the Monte Carlo work of Ostrowski 
1991). We find that in our application to non-relativistic shocks, the choice of scattering mode 
is largely immaterial to the resultant distributions; we expect this not to be so for relativistic 
shocks where the modes generate vastly different particle anisotropics. We also compare non- 
linear (Section 4.2) results with test-particle ones where the non-thermal particles do not modify 
a discontinuous shock (Section 4.1), finding that, as in our earlier work on plane-parallel shocks, 
some spectral curvature arises in high Mach number shocks where a large fraction of the partial 
pressure resides in the non-thermal population. An outline of the Monte Carlo method is given in 
Section 2, followed in Section 3 by flux conservation considerations and the associated scheme for 
iterative determination of the modified shock flow and field profiles. The spectral and flux results 
comprise Section 4, culminating in a presentation of acceleration efficiencies and discussion of the 
results. 



2. THE MONTE CARLO METHOD 

The Monte Carlo technique for describing particle acceleration at plane shocks has been de- 
scribed in previous papers (e.g. Ellison, Jones and Reynolds 1990; Jones and Ellison 1991; Baring, 
Ellison and Jones 1993; Ellison, Baring and Jones 1995), and is essentially a kinematic model 
that closely follows Bell's (1978) approach to diffusive acceleration. The simulation follows indi- 
vidual particles as they traverse a background "plasma" consisting of an average bulk flow and 
magnetic field; the flow velocity and magnetic field consist of a grid of values from far upstream 
to far downstream with a subshock positioned at x = . Our subshock is a sharp discontinuous 
substructure of the overall shock, defining the conventional boundary between (infinite) upstream 
and downstream regions. Strictly it should be no sharper than the smallest diffusion scale (i.e. the 
gyroradius of thermal particles), however we require it to be abrupt for the purposes of expedience 
in the simulation. Particles are injected far upstream of the shock with a thermal distribution at 
temperature T\ , mimicking for example solar wind ions (as in applications to the Earth's bow 
shock), and are allowed to convect in the flow and scatter, crossing the shock a few or many times 
before they eventually leave the system either far downstream or beyond an upstream free escape 
boundary (FEB). They are moved one at a time according to a prescribed scattering law, defined 
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below, in a test-particle fashion until they exit the simulation. In cases where efficient acceleration 
arises, feedback of the accelerated population leads to significant smoothing of the shock profile 
and heating of particles occurs in the foreshock region; discussion of such non-linear aspects of the 
simulation are deferred to Section 3 below. 

Following our previous treatments of oblique shocks (Baring, Ellison and Jones 1993; Ellison, 
Baring and Jones 1995) particle convection is performed in the de Hoffmann- Teller (HT) frame (de 
Hoffmann and Teller 1950), a frame where the shock is stationary, the fluid flow u is everywhere 
parallel to the local field B , and the electric field is u x B = everywhere. The HT frame of 
reference is therefore particularly convenient because of the associated absence of drift electric fields: 
particle trajectories are then simple gyrohelices and the description of convection is elementary. 
Furthermore, it follows from the mere existence of an HT frame that the so-called shock drift 
mechanism is inseparable from, and intrinsically part of, the Fermi acceleration process (e.g. Drury 
1983; Jones and Ellison, 1991) and is therefore automatically included in our Monte Carlo technique 
since particle motion is followed in the HT frame. 

While field and flow directions in the HT frame are uniquely defined downstream of the shock, 
the non-linear nature of this work yields a spatial variation of u and B , due to the compressive 
effects of the accelerated population. This variation is accommodated using a grid-zone structure 
that was implemented in many earlier versions of the Monte Carlo technique; each zone contains 
uniform field and flow, with discontinuities at the boundaries satisfying the Rankine-Hugoniot 
conditions discussed in Section 3 below. The grid-zone boundaries therefore mimic mini-shocks, 
with the subshock defining a particular grid boundary; a depiction of this grid structure is given 
in Baring, Ellison and Jones (1992). The spatial resolution of the grid can be adapted at will, but 
in our applications, we require it to be finer at some distance upstream than the typical mean free 
path of particles that penetrate to that distance from the shock. Note that the ability to define 
(for example via the Rankine-Hugoniot conditions) an HT frame with u x B = on both sides of 
a grid point implies, by spatial extension, that the HT frame is uniquely defined throughout the 
flow, regardless of flow and field compression and deflection upstream. In addition, note that even 
though the u x B electric field is transformed away by going to the HT frame, charge separation 
electric fields are not and these have not been included in our model, being beyond the scope of 
the present work. 

While particle transport is monitored in the HT frame, all measured quantities such as particle 
distributions and momentum and energy fluxes are output in the normal incidence frame (NIF), 
which is the frame where the shock is also stationary, but is defined such that the flow far upstream 
(i.e. where it is uniform to infinity) is normal to the shock plane. A simple velocity boost with a 
speed of v HT = u\ tan ©Bni parallel to the shock front effects transformation between the NIF and 
HT frames, where 0Bni is the far upstream angle the magnetic field makes with the shock normal 
and Mi is the far upstream flow speed in the frame where the shock is at rest. Note that hereafter, 
the index '1' will indicate far upstream values and the index '2' will indicate far downstream values 
well away from the smooth shock transition. A depiction of the NIF geometry is given in Figure |l] 
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for the specific case of unmodified shocks; modified shock geometry extends this to include piecewise 
increments of u and B . In the normal incidence frame, the — x-axis defines the shock normal, and 
the senses of the other axes are as in Figure [l|. The results of this paper are restricted to highly 
subluminal cases where v UT <C c , the speed of light, since our Monte Carlo technique has not yet 
been generalized to include relativistic effects in oblique shocks. 

The simulation is an orbit code, where particle propagation is performed by following the 
gyromotions exactly, as in the test-particle work of Ellison, Baring and Jones (1995) and a diversity 
of works in the literature (e.g. Decker 1988; Begelman and Kirk 1990; Ostrowski 1991; Takahara 
and Terasawa 1991). The position of particles is incrementally updated on a timescale 5t , which 
is a small fraction of a gyroperiod, i.e., 5t <C r g = mc/ (QeB) , where m and Q are the particle's 
mass and charge number respectively, and e is the electronic charge. A particle in a particular 
grid-zone is moved in a helical orbit determined by the magnetic field and bulk flow velocity for 
that grid position. 



2.1. Particle Scattering 

Having outlined the procedure for convection, here we describe our prescription for particle 
scattering, which is somewhat more involved. After each time step 5t a determination of whether 
the particle should "scatter" or not is made using a scattering probability P S cat = St/t c where the 
collision time, t c , is given by X/v , A is the mean free path, and u < c is the particle speed, both 
measured in the local fluid frame. This prescription yields an exponential pathlength distribution. 
The scatterings, presumably off magnetic irregularities in the flow, are assumed to be elastic in the 
local plasma frame so that monotonic energy gains naturally arise as a particle diffuses back and 
forth across the shock due to the converging nature of the flow. It proves convenient to scale the 
mean free path by the gyroradius, introducing a model parameter r] that is the ratio of the two 
quantities (following Jokipii 1987; Ellison, Baring and Jones 1995): 

A = 7/r g or K|| = ^r]r g v , (1) 

where km is the diffusion coefficient parallel to the ambient local field. It follows that the collision 
time satisfies 

A rjr g mc , n s 

v v ' QeB y ' 

Generally, rj is a function of energy, however in this paper it is assumed to be a constant inde- 
pendent of position and energy (also following Jokipii 1987). For this choice, the collision time is 
independent of particle energy, a convenient simplification that can be easily generalized to some 
other dependence of A on r g (see Ellison, Mobius, and Paschmann 1990). Note also that A is 
hence implicitly inversely proportional to B . 

The assumed constancy of rj is the most important approximation we make, since all of the 
complicated plasma physics of wave-particle interactions is incorporated in equation (|l|). Conve- 
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nience aside, there are sound reasons for choosing this simple relation. First of all, as long as 
electrostatic effects are neglected (they are omitted from our treatment), the gyroradius is the 
fundamental scale length of a particle at a particular energy and the mean free path can be ex- 
pected to be some function of this parameter. Second, if the plasma is strongly turbulent with 
5B/B ~ 1 , as is generally observed in space plasmas, the large-scale structures in the magnetic 
field will mirror particles effectively on gyroradii scales (i.e. for rj rs 1 , the Bohm diffusion limit; 
Zachary 1987). Third, and most important, spacecraft observations suggest that A « r g in the 
self-generated turbulence near the earth's bow shock (e.g. Ellison, Mobius, and Paschmann 1990). 
Fourth, hybrid plasma simulation results also suggest that the mean free path is a moderately 
increasing function of r g (Giacalone, Burgess, and Schwartz 1992; Ellison et al. 1993). Note that 
while the exact form for A will surely have quantitative effects on the injection rate and all other 
shock characteristics, the most important qualitative effects should be well modeled as long as a 
strongly energy dependent diffusion coefficient is used. Employing a realistic Kit with a strong 
energy dependence is essential because of the intrinsic efficiency of shock acceleration. If ku is 
indeed energy dependent, the highest energy particles with large fractions of energy and pressure 
have very different scales from thermal particles, leading to the spectral curvature that appears in 
the simulated distributions (see Section 4). 

The simulation employs two complementary types of scattering modes, namely large-angle 
scattering and pitch-angle diffusion. For each mode, elasticity of scattering is imposed, which 
amounts to neglecting any recoil effects of wave production on the particles and hence that the 
background scattering centers (i.e. magnetic irregularities) are frozen in the plasma. This approx- 
imation is generally quite appropriate but becomes less accurate for low Afvenic Mach numbers 
when the flow speed does not far exceed the Alfven wave speed. In large-angle scattering (LAS), 
a particle's direction is randomized in a single scattering event (on a timescale of t c ) and the new 
direction is made isotropic in the local plasma frame. Such quasi-isotropic scattering is adopted 
in most of our earlier simulation work (e.g. see Jones and Ellison 1991; Baring, Ellison and Jones 
1993; Ellison Baring and Jones 1995), and is intended to mimic the effect of large amplitude field 
turbulence on particle motions. Such turbulence is present in both plasma simulations (e.g. Quest 
1988; Burgess 1989; Winske et al. 1990) and observations of shocks in the heliosphere (e.g., Hoppe 
et al. 1981). 

The second scattering mode we employ is pitch-angle diffusion (PAD), as used in the diverse 
works of Decker and Vlahos (1985), Decker (1988), Kirk and Schneider (1987) and Ostrowski (1988, 
1991). In this mode, the direction of the velocity vector v is changed by a small amount after 
each time step, St , rendering the scattering process more "continuous" than large angle collisions, 
and more appropriate to physical systems with small levels of field turbulence. Here we adopt the 
procedure detailed in Ellison, Jones, and Reynolds (1990) for determining the maximum amount 
the pitch angle 9 B = arccos{v-B/|v| |B|} can change after each 5t ; our procedure is summarized in 
the Appendix. A comparison of the simulation results for these two modes is one goal of this paper, 
motivated by an expectation that they could, in principal, produce different injection efficiencies at 
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modified shocks. This expectation is partly based on the spectral differences observed between LAS 
and PAD applications to unmodified relativistic shocks (e.g. Ellison, Jones, and Reynolds 1990), 
where the two modes generate significantly different particle anisotropics; distribution anisotropies 
are indeed relevant to the injection problem considered here. 

As a particle convects, the simulation tracks both its position and the position of its gyrocenter. 
After a scattering occurs and a new direction is obtained for its velocity vector, a new gyrocenter 
is calculated. This shift of the gyrocenter means the particle is now gyrating around a different 
field line and diffusion across the field has occurred; the new field line is within 2r g of the one the 
particle was circling before the scattering. Such cross field diffusion is an integral part of diffusive 
acceleration at oblique shocks (e.g. Jokipii 1987, Ellison, Baring, and Jones 1995), and its presence 
is required in the Monte Carlo simulation in order to match spacecraft observations of particle 
spectra associated with interplanetary shocks (Baring et al. 1995). Ellison, Baring and Jones 
(1995) showed that this scheme for cross-field diffusion together with the assumption contained in 
equation (|l]) is equivalent to a kinetic theory description of diffusion (e.g. Axford 1965; Forman, 
Jokipii and Owens 1974), where the diffusion coefficients perpendicular to (k±) and parallel to 
(/cm ) the field are related via k± = /cm/(1 + rj 2 ) . The parameter r\ in equation (|l]) then clearly 
determines the strength of the scattering and when 77 ~ 1 , kj_ ~ km , the so-called Bohm limit 
where particles diffuse across the magnetic field as quickly as they move along it. The properties 
of highly oblique and quasi-parallel shocks tend to merge when the scattering is strong. 



2.2. Grid Zone, Free Escape and Downstream Return Boundaries 

When a particle crosses a grid zone boundary, the values of the bulk flow velocity and the 
magnetic field (both magnitude and direction for these vector quantities) change and a new gyro- 
radius, gyrocenter, and phase in the gyro-orbit are determined, as outlined in the Appendix. The 
particle then acquires a new gyromotion with subsequent convection along the new field direction. 
At a crossing of a grid zone boundary, we adopt the standard requirement (e.g. see Terasawa 
1979; Decker 1988; Begelman and Kirk 1990) in orbit calculations that the momentum vector of 
the particle is conserved in the de Hoffmann- Teller frame, since there is no electric field in the 
shock layer in this frame. While this implies conservation of energy in the HT frame (i.e. between 
scatterings), particle energies do change at the discontinuity in the NIF (e.g. Toptyghin 1980) 
due to the presence of drift electric fields, a manifestation of the transformation between frames. 
This transmission criterion differs from the imposition of magnetic moment conservation that was 
made in earlier applications of our Monte Carlo technique (e.g. Baring, Ellison and Jones 1993); 
in the present paper we make no assumption concerning the magnetic moment of a particle at a 
grid point or anywhere else. Particles may of course 'reflect' at any zone boundary or be trans- 
mitted depending on their phase and pitch angle after a number of gyrations in the vicinity of the 
boundary. 

There are two limiting boundaries to the simulation region, the free escape boundary (FEB) 
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and the downstream return boundary. Our introduction of an upstream FEB facilitates modelling 
of the finite extent, for example through geometrical curvature, of real shocks. Escape is naturally 
expected in real systems since the region of influence of a shock on its environs is finite, and the 
level of shock-generated turbulence diminishes to background levels at sufficient distances from the 
shock. An escape boundary is most relevant to the upstream region (i) because the direction of 
convection renders the downstream region more spatially uniform, and (ii) because the upstream 
region is usually on the convex side of the shock (e.g. supernova remnant shocks, the Earth's bow 
shock), although particles can escape from the downstream side as well. The inclusion of such a 
free escape boundary is also motivated on theoretical grounds. The fundamental point is that for 
Fermi acceleration in steady-state shocks, escape must occur for Mach numbers above some critical 
value. This has been fully documented (e.g. Eichler 1984, 1985; Ellison and Eichler 1984; Jones and 
Ellison 1991), where it is observed that within the context of the non-linear acceleration model, 
the Fermi acceleration/hydrodynamics coupling becomes unstable for shocks of Mach numbers 
above a few, and leads to singularities in the energy density when particle escape is suppressed. 
Finite solutions are achievable when an upstream FEB is introduced, since its presence causes the 
acceleration process to truncate at the highest energies. In the case where the diffusion coefficient 
increases with energy, the FEB produces a distribution which falls off approximately exponentially 
at an energy where the upstream diffusion length is on the order of the distance from the FEB to 
the shock. In the simulation, particles are removed just before they scatter for the first time on the 
upstream side of the FEB; this choice leads to a spatial smearing of the effects of the FEB on the 
scale length of the mean free path of the escaping particles, i.e. on length scales comparable to the 
distance d FEB between the FEB and the shock. In the results presented here, the FEB is chosen 
close enough to the shock to guarantee that all particles in the simulation remain non-relativistic; 
the domain of acceleration to relativistic energies and also relativistic shock scenarios are deferred 
to future work. The dynamical consequences of the FEB are discussed in detail in Section 4 below. 

While the upstream region in the simulation is finite, delimited with a FEB, we model an 
infinite downstream region with a probability of return calculation beyond a downstream return 
boundary (DRB); this spatial border renders the simulation finite in time. Beyond the DRB, which 
is maintained more than a scattering length downstream of the shock, the spatial diffusion properties 
are treated using appropriate statistical probabilities. If the position of a particle (as opposed to its 
guiding center) is followed, then the probability of return to the upstream side of the DRB assumes 
a simple form. If the flow is uniform with a component of velocity u x 2 perpendicular to the DRB 
(in our case, this direction is also perpendicular to the shock plane) and particles of speed v F in 
the frame of the flowing plasma are also isotropic in that fluid frame, then the probability, P re t , 
that a particle which crosses some arbitrary y-z plane will return to the upstream side of that 
plane, is 

Piet= (?hL^)\ (3 ) 
\v F + u x2 J v ; 

While this calculation has been done many times (e.g. Bell 1978; Drury 1983; Jones and Ellison 
1991; and Ellison, Baring, and Jones 1995) we emphasize that equation (||) is fully relativistic 
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(Peacock 1981) and holds regardless of the orientation of the magnetic field or the flow. The 
principal requirement for the validity of equation (||) is that the particles are isotropic in the local 
fluid frame, a condition that is satisfied since the DRB is at least a scattering length downstream 
of the shock. Hence while equation (|3|) can be used downstream where the flow is uniform for any 
Vp ^ u x 2 , it cannot be used at the shock where the flow speed changes unless u x 2 <C v F . The 
particle speed, v F , must also remain constant during the time a particle spends downstream from 
the y- z plane, a natural consequence of our elastic scattering assumption. The decision of return 
(or otherwise) is made via a random number generator. Particles which do return, must be injected 
back across the y-z plane with properly flux-weighted ^-components of velocity, pitch angles, and 
phases. The determination of these, along with a detailed derivation of equation (||), are given in 
the Appendix. Note that the DRB is not only a feature of our Monte Carlo simulation, but is also 
used in hybrid plasma simulations of shocks (Bennett and Ellison 1995). 

3. FLUX CONSERVATION RELATIONS AND SHOCK MODIFICATION 

Before presenting the results of our modified shock simulations, it is instructive to review 
the elements of non-linear shock hydrodynamics and our procedure for determining the fluid flow 
and magnetic field spatial profiles that simultaneously conserve all relevant fluxes and are also 
self-consistent products of the Fermi acceleration mechanism. 

3.1. Flux Conservation Relations 

The starting point for these considerations is the well-known one-dimensional, steady-state, 
magnetohydrodynamic conservation relations (i.e., the Rankine-Hugoniot (R-H) jump conditions) 
for an infinite, plane shock lying in the y-z plane (see the geometry in Figure |l]). Variations of all 
quantities occur only in the x -direction and these equations are written in the normal-incidence- 
frame [NIF]. The notation is that of Decker (1988), with the square brackets representing differences 
between quantities far upstream (with the '1' subscript) and downstream (with the '2' subscript) of 
the shock, however the origin of the forms used here is based on the presentation on p. 56 of Boyd 
and Sanderson (1969). For a magnetic field strength of B , if u is the bulk speed of the plasma, 
the purely electromagnetic equations (i.e. Maxwell's equations) are 

ml = o , (4) 

which defines a divergenceless magnetic field (remember that our system has dB y /dy = = 
dBz/dz ), and 

[u z B* - UxB z ]l = , (5) 

which expresses (since c V x E = — <9B/<9i = ) the uniformity of the tangential electric field across 
the shock. The hydrodynamic equations are as follows: the mass flux equation corresponding to 
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the x -direction is 



[pu x ]i 







(6) 



where p is the mass density; the equations for the flux in the x -direction of the x and z compo- 
nents of momentum are 
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respectively, where P xx and P xz are the appropriate components of the pressure tensor. Finally, 
the energy flux in the a; -direction satisfies 
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Here 7 is the ratio of specific heats, which enters via the thermal contribution to the energy density 
using the equation of state; we set this equal to 5/3 in this paper, since only non-relativistic particles 
appear in the simulation results presented. Equations (^)-(^) neglect so-called gradient terms that 
are spatial diffusion contributions that arise from non-uniformity of the fluid flow and magnetic 
field profiles. Note also that equations (0), ®, and (g) approximate the respective parallel shock 
relations for high Alfvenic Mach numbers (i.e. where the field is dynamically unimportant), whereas 
equation (||) remains important regardless of the Mach number. 

In equation @ we have added the term, Q csc , to model the escape of particles at an upstream 
free escape boundary (FEB). As mentioned above, a FEB causes the acceleration process to truncate 
as particles leave the system producing important dynamical effects since the escaping energy, and 
therefore pressure, results in an increase in the compression ratio of the shock (see Ellison, Mobius, 
and Paschmann 1990 for a discussion of the effects of such a term in the R-H relations). The 
escaping energy flux, Q esc , is taken to be constant for the far downstream region and zero for 
the region far upstream (i.e. several mean free paths) of the FEB, and varies most rapidly in the 
neighbourhood of the FEB. We assume that the concurrent escaping momentum and mass fluxes 
are small and neglect them in equations @, @), and This is a good approximation if the 
particles that escape have speeds such that v csc S> u x \ (see Ellison 1985), a situation that is always 
realized in the simulation results presented here. 

The appearance of different components of the pressure tensor in the Rankine-Hugoniot rela- 
tions is requisite for the Monte Carlo simulation since we do not assume that particles are isotropic 
in any frame. Normally, implementations of the conservation equations in astrophysical or helio- 
spheric applications (e.g. Decker 1988) are restricted to scenarios where the plasma is isotropic in 
the local fluid frame, in which case = P yy = P zz = P and the off-diagonal terms of the pressure 
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tensor are zero. However, our system generates anisotropic plasma in all frames of reference, due 
to the non-uniformity of the flow combined with the self-consistently determined Fermi acceleration 
of the particles. In this paper, for the sake of simplicity, we assume that the plasma is gyrotropic 
but anisotropic in the local fluid frame in the flux equations. 

When an isotropic population of particles is convected across a velocity discontinuity and 
then subjected to isotropic scattering, compression of the plasma is close to, but not perfectly, 
gyrotropic: asymmetric phase sampling at the discontinuity yields non-uniform phase distributions 
prior to scattering. Plasma isotropy in the new local fluid frame is attained only after many 
scattering lengths; in fact adjustment to true isotropy is never achieved in our non-linear Monte 
Carlo treatment because the scale length of velocity (and directional) changes in the fluid flow is 
always comparable to the mean-free-path of the particles comprising the flow. Gyrotropy is a good 
approximation that is also expedient because it yields a diagonal pressure tensor Pj in the de 
Hoffman- Teller and fluid frames for co-ordinate systems with one axis aligned along the magnetic 
field; since the diagonal components generally differ, oblique shocks lead to non-zero off-diagonal 
components: P xz = Pzx 7^ . A discussion of the generation of such terms is presented in the 
Appendix, specifically focusing on the details of the derivation of the pressure terms of the energy 
flux in equation @; there the coefficient of P xz is alternatively expressed in terms of pressure 
components parallel to and orthogonal to the local field. Note that non-zero P X y — ^yx also arise 
in hydrogenic plasma flows in conjunction with out-of-the-plane components of the electric field 
(Jones and Ellison 1987, 1991); such off-diagonal terms apply only to quantities in the y-direction, 
and therefore are irrelevant to the considerations of this paper. 



The flux equations are used in the simulation in dimensionless form, scaling by relevant 
upstream quantities. In the NIF, the far upstream flow is taken along the shock normal, i.e. 
u z \ = . We defi ne the flow velocity u\ = u\ x , the magnitude of the far upstream mag- 
netic field B\ = \J Bf x + P>1 z , and a far upstream plasma density p\ . These specify a far 
upstream Alfvenic Mach number: M A \ = u\/v A \ or M\\ = (4-7T p\u\) / B\ , where the Alfven 
speed is v A \ = Bi/(4npi) 1 ^ 2 , and a far upstream sonic Mach number, M| x = p\u\/ {^P\) , where 
Pi = nik B Ti is the far upstream (isotropic) pressure. Here n\ and T\ are the number density and 
temperature far upstream, and k B is Boltzmann's constant. These upstream parameters, along 
with @Bn ; define the key input for the simulation runs. Using these definitions, one can write 
equations (||) and (||)-(||) in a dimensionless form at any position x : 



3.2. Flux Scalings and Formalism 




(10) 



-13- 



defines the uniformity of tangential electric field, 
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(12) 



u' x (x)B' z (x)-u' z (x)B' xl +Q' esc = F' enl , (13) 



rearranges the energy flux equation. All primed quantities are dimensionless, using the notation 
v! = uju\ , B = BjB\ , P' = P/{piu\) and Q esc = Qesc/ {p\u\) ■ Note that 7 = 5/3 is a constant 
throughout the flow for the non-relativistic applications here. The constancy of the x-component 
of magnetic field has been used to substitute B x (x) = B xl , and we have also used mass flux 
conservation 

p(x)u x (x) = p\U\ = constant , (14) 

in these equations. As mentioned above, the escaping momentum and mass fluxes are of progres- 
sively smaller orders in u x i/v esc than Q' esc , and therefore are neglected. 



The far upstream fluxes on the right hand sides of equations 
by the input shock parameters: 
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(15) 



Note that in all of the examples described below, we have set the electron temperature equal to 
zero. A finite electron temperature can be modeled with our procedure and is necessary when 
fitting spacecraft data (as done in Ellison, Mobius, and Paschmann 1990), but is not necessary 
for the discussion given here. We assume, as is generally done, that the ions dominate the shock 
structure and the addition of electrons has little effect other than changing the Mach number. 

If Q csc and P xz are both assumed to be zero at all x , the four unknowns in equations (|l0|) -(p^); 
u x (x) , u' z (x) , P4x( x ) > an d B z (x) , can be obtained at every position x . This is just the standard 
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situation of a discontinuous shock, and these Rankine-Hugoniot relations are analytically solvable 
(e.g. see Decker 1988) for the shock compression ratio, r = u\/u x 2 . However, in the modified 
collisionless shocks considered here, the non-thermal component of the particle distribution that 
is generated by particles crossing the shock more than once contributes significantly to the total 
pressure of the system, and Q esc will not, in general, be zero. In fact, Q esc will have different 
values at various locations, and cannot be determined before the shock structure is known, making 
a direct solution of equations (|i~0|)-(13) impossible; this is the inherent non-linearity in the problem 



even if isotropy is assumed, defined by the coupling between the acceleration process and the flow 
hydrodynamics. 

In our approach, we iterate to achieve a solution for the velocity and field shock profiles by 
varying v! x (x) , u' z (x) , B z (x) and the overall compression ratio, for successive simulation runs 
(each accelerating particles and generating non-thermal distributions) until equations (|i~0|)-(|i"^) 
are satisfied at every x . The overall compression ratio depends on Q esc far downstream from 
the shock and is determined by our solution. When this value is consistent with equation (|i~3|), a 
complete solution to the non-linear acceleration problem is obtained, satisfying equations (||)-(||) 
at all positions. The details of the iterative procedure follow. 



3.3. Iteration of the Shock Profile 

The iteration of the shock profile is done in two stages. We first choose the overall compression 
ratio (normally the R-H value for the first iteration) and, using this ratio, we iterate the shape of 
the profile. As individual particles move through the shock, the momentum and energy fluxes are 
calculated at each grid zone boundary. We therefore obtain the quantity 

FUx) = u' x (x) + PUx) + |^ (16) 

at each boundary, where u' x (x) and B z (x) are the current values for the shock structure. From 
this we compute the pressure P^^x) and calculate a new x -component of the flow speed from 



equation (|ll|), 



u?{x) = F'^ - PUx) - fj^ (17) 

such that the momentum flux will equal the constant far upstream value, F^ , at all x . Replacing 
v! x (x) with u x (x) , we solve equations (|i~0|) and ( |l2|) for new values of u z (x) and B z {x) . To speed 
convergence, before running the next iteration we smooth this new profile, force it to be monotonic, 
average it with the previous profile, and scale the profile by setting all downstream values of the 
x -component of flow to ii x 2 and the far upstream value to u\ x . Since no wave physics is employed 
in our description, we do not attempt to model anything other than a monotonic decrease in 
flow speed and a monotonic increase in the magnitude and obliquity of B from upstream to 
downstream. Alternatively, we can calculate both P^ix) and P xz (x) in the simulation and obtain 
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the new prediction for u^ N (x) from equation (|l3|), but in either case, our procedure rapidly and 
stably converges. 

With this new shock profile we repeat the simulation by again injecting particles far upstream 
from the shock and propagating them until they leave at either the FEB or the probability of 
return plane. Our algorithm converges rapidly (within a few steps: see the examples in Section 4 
below) to values of u' x (x) , u' z (x) , and B z (x) , which no longer change significantly with subsequent 
iterations. However, in general, this profile will not simultaneously conserve momentum and energy 
fluxes unless a compression ratio consistent with the escaping energy flux Q eS c has been chosen. 
The second stage of the iteration process is to successively choose new overall compression ratios 
(larger than the R-H value), each time repeating the iteration of the profile shape, and continue 
until the momentum fluxes and the energy flux (with Q esc added) are constant everywhere. Thus, 
within statistical limits, a shock profile and overall compression ratio that are consistent with the 
Fermi acceleration process are obtained. 

4. RESULTS 

Oblique shocks are highly complex, even in the steady state and in plane geometry, and several 
parameters control the dissipative processes as well as the injection from thermal energies into the 
Fermi acceleration mechanism. These far upstream parameters include the magnetic field strength, 
B\ , the obliquity, ©Bni > the temperature, T\ , the number density, n\ , and the shock speed, 
u± , all of which are determined by the ambient upstream conditions and can, in principle, be 
determined by observations of a given physical system. The size of the acceleration region is also 
an observable (for example the radius of a supernova remnant shock), and we model it using the 
distance d FEB between the upstream free escape boundary and the shock. However, the 'size' of 
the shock in units of mean free paths is very important and this will depend on the scattering law 
we assume. This requires the introduction of another parameter, r\ , the ratio of the mean free path 
to the gyroradius, via equation (||). The value of t] , which determines the amount of cross- field 
diffusion, depends on the highly complex plasma interactions that occur in the shock environs; the 
prescription in equation (||) is a simple but insightful way to model these plasma processes. 

Another "variable" results from the inclusion of two extreme modes of scattering, namely large- 
angle scattering (LAS) and pitch-angle diffusion (PAD). While more complicated scattering models 
can be used, we believe these contain the essential physics of plane shocks and yield important 
information on the nonlinear processes linking shock structure and particle acceleration. The type 
of scattering we employ and r\ are free parameters and cannot be determined in our model except by 
comparison with observations of space plasma shocks or 3-D plasma simulations. Hence (replacing 
B\ and T\ with M A \ and M s \ ), there are seven parameters in our model: Mai > M s i , n\ , 
©Bni , ui , d FEB , and r\ , together with the choice of the type of scattering (either LAS or PAD). 
In all of the following examples we use a shock speed of u\ = 500 km sec -1 and a far upstream 
number density of n\ = lcm~ 3 , which define physical scales for our system that are more or less 
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appropriate for astrophysical shocks. The spatial scales of the results presented here are all in units 
of the "low-energy mean free path" Ao , which is defined as the mean free path [see equation (|l])] of 
a proton with speed u% in the upstream magnetic field, i.e., Ao = i]m^u\cj {eB{) , where rj (> 1) 
remains an adjustable parameter for each simulation run. 



4.1. Test-Particle Examples 

The simplest acceleration results obtainable from the Monte Carlo simulation are for test- 
particle cases where the shock profile is uniform on either side of the subshock; this limit corresponds 
to the first run in the iteration sequence described in Sections 3.2 and 3.3, and has been studied in 
detail in Ellison, Baring and Jones (1995). Several interacting elements of the code, including shock 
(or grid zone) crossings and the probability of return calculation (which includes the flux- weighting 
of momenta of the returning particles: see Section A. 3 in the Appendix) must be implemented 
properly in order to yield the well-known test-particle acceleration power-law. The Fermi power- 
law is achieved when particle speeds v far exceed the HT flow speed ii x i/ cos Beni • From analytic 
calculations (e.g. Drury 1983; Jones and Ellison 1991), the spectral index a of the power-law then 
depends only on the compression ratio r = u\lu<z^ regardless of the obliquity or other plasma 
parameters. The compression ratio is determined from equations ((§H@ with Q esc = -P X z = 
(see also Decker 1988). For non-relativistic particle energies and shock speeds the test-particle 
distribution is: 

^cxE^ , a = , r + 2 , , (18) 
dE 2(r - 1) y ' 

where dJ/dE is the number of particles in units of (cm 2 -sec-steradian-keV)~ 1 , i.e. is an omni- 
directional flux (see Jones and Ellison 1991). The reproducibility of this form is a powerful tool for 
debugging the portions of the code that are directly related to the transport and acceleration of 
particles. It is instructive to review test-particle distributions before proceeding to our results for 
the non-linear problem. 

In Figure ^ we show spectra calculated with a discontinuous shock for three different sets of 
parameters and for both large-angle scattering (solid lines) and pitch-angle diffusion (dotted lines). 
Note that the examples labelled (c) are multiplied by 0.01 for clarity of display. All spectra here and 
elsewhere are omni-directional, calculated several Ao (i.e. "thermal" mean free paths) downstream 
from the shock in the normal incidence frame, and are normalized to one particle per square cm 
per second injected far upstream. In all cases we have used u\ = 500km sec -1 , and rii = 1cm -3 ; 
the (sonic and Alfvenic) Mach number 20 cases here use B\ = 1.15 x 10 -5 G and T\ = 4.54 x 10 
K, while the Mach number 3 case uses B 1 = 7.64 x 10 -5 G and T x = 2.02 x 10 4 K. The other 
parameters are given in the figure. No free escape boundary is included in these test-particle cases, 
being necessary only for the non-linear acceleration problem. 

Note that the compression ratio varies slightly between the two high Mach number examples, 
being r = 3.96 for the Beni = 30° case and r = 3.94 for the 60° case; this occurs because of 
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a slight dependence of r on obliquity for non-infinite Mach numbers. For the low Mach number 
example (c), the shock is quite weak with r = 2.5. The Fermi power-laws obtained from equa- 



simulation does reproduce the Fermi power-law index at high energies for a wide range of shock 
parameters. Comparison of examples (a) and (b) with similar compression ratios but quite different 
0Bni j supports the fact that the Fermi spectral index a is determined solely by r . 

The next most striking feature of these plots is that the two modes of scattering produce very 
little difference in the spectra. This difference is largest for portions of the ©Bni = 60° spectrum 
(b) between thermal energies and about 100 keV. There the distribution for the LAS mode [the 
solid line in (b)] is somewhat noisy due to comparatively poor statistics in high rj , high obliquity 
runs (i.e. weak injection: see Ellison, Baring and Jones 1995), and this noise may obscure some 
underlying structure that can arise from large energy boosts in individual shock crossings. At 
these energies, the particle speed does not far exceed the HT frame flow speed u x i/cos 0Bnl > 
leading to significant anisotropies in the particle population, and more importantly, measurable 
differences in the degree of anisotropy produced in PAD and LAS modes. Since such differences in 
angular distributions are responsible for observed differences between LAS and PAD applications 
to unmodified relativistic shocks (e.g. see Ellison, Jones, and Reynolds 1990 for a comparison of 
the modes, and Kirk and Schneider 1987; Ostrowski 1991 for PAD cases), it is not surprising that 
spectral differences should appear here at suprathermal energies for highly oblique shocks. Clearly, 
Figure || shows that for this intermediate obliquity case, the mode of scattering has virtually no 
effect on the resultant spectrum at either thermal or the highest energies, and therefore that the 
scattering mode plays little role here in determining the efficiency of acceleration, i.e. the fractional 
energy deposited in high energy particles. This is not surprising, because of the low value of rj 
here: we naturally expect that strong scattering (near the Bohm diffusion limit) will destroy any 
sensitivity of the acceleration process to the shock obliquity or distribution anisotropies, and hence 
the type of scattering. Later, we shall see that larger (but still relatively small) differences between 
the scattering modes occur for weak scattering (i.e., large rj). 

Also evident in Figure || is the strong effect the input parameters have on injection efficiency: 
the @Bn = 60° spectra fall an order of magnitude below the @Bn = 30° spectra at high energies 
even though they both obtain the same power law index. Increasing either ©Bni or rj will result 
in decreased injection efficiency. These effects were detailed in Ellison, Baring, and Jones (1995), 
where an anti-correlation between acceleration time and efficiency of acceleration in test-particle 
shocks was observed. We note that existing analytic predictions for the transition between the 
thermal peak and the high energy power-law, i.e., the injection efficiency, require ad hoc parameters 
additional to and independent of those made for the shock structure to connect the thermal gas to 
the cosmic ray population (e.g. Zank, Webb, and Donohue 1993; Kang and Jones 1995; Malkov and 
Volk 1995). The advantage of our model is that the single relation [equation (yj)] controls the shock 
structure, the absolute injection efficiency, and, in fact, the entire shock solution. This makes it 
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straightforward to compare model predictions to observations and to infer plasma properties (such 
as the level of turbulence, the correctness of the elastic scattering assumption, etc.) from these 
comparisons, an attractive feature. Properties of upstream particle distributions are deferred to 
the discussion of non-linear results in the next subsection, though test-particle spectra upstream of 
oblique shocks were presented in Baring, Ellison and Jones (1994). 

4.2. Examples Showing Iteration of Shock Profile 

For our next examples, we compute the self-consistent smooth shock profile beginning with a 
low Mach number case, i.e., M A \ = M s \ = 3 , Oeni = 30° , d FEB = — 50Ao , and r) = 2 , yielding a 
plasma /3 of (3i = 1.2 . This corresponds to a weak shock, typical of interplanetary shocks observed 
in the heliosphere (e.g. see Burton et al. 1992, Baring et al. 1995). To reiterate, in the results that 
follow, all lengths are measured in units of Aq , which is the mean free path rjr g i of a proton of 
gyroradius r g i = m p u\c/ {eB\) , i.e. with the speed u\ in the upstream magnetic field. 

In the left panels of Figure [3| we depict the average flow speed, u x , the flux F^(x) of the 
x-component of momentum, and the energy flux F en (x) , all normalized to far upstream values, 
for several iterations starting from a discontinuous shock (light solid line) and yielding the final 
profile (heavy solid line). These iterations were done using LAS and an overall compression ratio 
of r ~ 2.7 which was determined in previous iterations on r . The Rankine Hugoniot compression 
ratio with Q csc = is r = 2.67 which is equal (within errors) of our 2.7 value. The convergence 
is quite rapid and the heavy solid lines (fourth iteration) show no further statistically significant 
change with additional iterations. Except for a departure of about 5% near x = , the flux of 
the x-component of momentum (middle panels) is constant for all x after the final shock profile 
has been obtained. The flux F^ zl of the z-component of momentum is generally small and less 
interesting; its profile (and those of u z and B z ) is not displayed for reasons of brevity. For the 
discontinuous shock (light solid lines), the momentum flux was clearly not conserved and rose to 
~ 140% of the far upstream value downstream from the shock. 

The escaping energy flux at the FEB (which is at — 50 Ao and not shown in the figure), is less 
than 1% of the far upstream value and does not influence the overall compression ratio significantly. 
The check on the consistency of the final profile is that the momentum flux and the energy flux, 
including the escaping flux, must both be conserved. When this is achieved (i.e. corresponding 
to the solid lines), we have a unique, self-consistent solution. The ~ 10% discrepancy in the 
energy flux near x = is most likely the result of the strong gradients in the shock and/or 
agyrotropic pressure tensor terms in the Rankine-Hugoniot relations which we have neglected from 
our flux considerations. This discrepancy decreases rapidly with increasing Mach number, and so 
is greatest for this present case. 

Figure [3| also shows the shock structure obtained with pitch-angle diffusion PAD (right panels) , 
all other shock parameters being the same as the LAS case. Within statistics, the two scattering 
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modes give identical results, both for the shape of the profile and the overall compression ratio. 
This feature is not surprising given that the test-particle results of the previous simulation bore 
this similarity out. 

In Figure || we show the distribution functions generated by the smooth shocks of Figure |3|. 
The spectra are omni-directional, measured downstream from the shock, and calculated in the 
shock (i.e. NIF) frame. As with the shock structure, the distribution functions obtained with the 
two scattering modes are identical within statistics. This is to be expected because r] = 2 is close 
to the Bohm diffusion limit where, as mentioned above, isotropic diffusion will naturally obscure 
the differences between PAD and LAS. There is a large difference, however, between the smooth 
shock results and the test-particle, discontinuous shock (dotted line), which was obtained using the 
same compression ratio (r = 2.7) computed in the self-consistent solution of the non-linear LAS 
simulation. The light solid line is the Fermi power-law expected from r = 2.7; the test-particle 
spectrum attains this result before the falloff at ~ 100 keV produced by the FEB. The discontinuous 
shock produces more efficient acceleration at low energies than the smooth shock, which follows 
from the nature of the smoothed shock profile: low energy particles feel the effect of the subshock 
whereas only the high energy particles sample the full compression ratio of r = 2.7 . This same 
property is responsible for the upward curvature of the non-linear spectra, which never attain the 
Fermi power-law, but flatten towards it before falling off due to the FEB. 

We have also obtained solutions (not shown) using exactly the same parameters as above 
except with a FEB at d FEB = — 10Ao • The smaller shock system causes a cutoff at a lower energy 
than seen in Figure || and the shock structure is correspondingly on smaller length scales. However, 
the self-consistent compression ratio is still r ~ 2.7 , consistent with a Q eS c = , as expected, since 
these low Mach number shocks put a small fraction of the available energy into energetic particles 
regardless of the shock size. 

As a more extreme example, we show in Figure |5| a high Mach number shock ( M s i = M A i = 20 ) 
with a much larger shock size, i.e., the FEB is at <i FE B = — 200Ao ■ The other input parameters are: 
«i = 500 km sec" 1 , n x = 1cm" 3 , 9 B ni = 30° , rj = 2 , B t = 1.15 x 10~ 5 G, and T x = 4.54 x 10 4 K, 
yielding (5 = 1.2 . We have only used the LAS scattering mode since the PAD results are essentially 
identical. Note that the self-consistent compression ratio used here is r = 5 (well above the R-H 
value of r = 3.96 ) which has been determined with previous runs not shown. The distance scale 
in Figure |5| is logarithmic for x < — WXq and linear for x > — WXq . 

There are several important features of this shock solution. In the first iteration with no shock 
smoothing, the momentum and energy fluxes are wildly non-conserved with both of them obtaining 
downstream fluxes almost 15 times as large as the upstream values. Despite this, the subsequent 
iterations converge rapidly and by the fourth iteration the momentum flux is conserved everywhere 
to within 5% of the upstream value (the first, second, third, and fourth iterations are shown by 
light solid, dashed, dash-dot, and heavy solid lines, respectively; the flat dotted line indicates the 
far upstream value). The effects from the anisotropic terms in the momentum and energy fluxes 
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are less noticeable here than in the previous low Mach number examples. As with the previous 
examples, the shock is smoothed out to the FEB, however the subshock here is considerably more 
distinct showing a sharp discontinuity between the flow just upstream from the subshock and the 
downstream flow. The width of the subshock is well within one Ao • 

As for the energy flux, it falls about 20% below the far upstream value due to the particles 
lost at the FEB. The escaping energy flux, Q esc , which is zero far upstream, falls rapidly around 
d FEB and then becomes approximately constant into the downstream region. This results in a 
compression ratio, obtained by iteration in previous runs, of r ~ 5 , compared to the R-H value of 
r = 3.96. A compression ratio of r = 5 implies Q' csc = 0.084 and a Q' csc /F^ nl ~ 0.17, and when 
this is added to the energy flux shown in the bottom panel of Figure [|, we have a self-consistent 
solution with all fluxes constant at all x to within a few percent. Larger escaping fluxes are 
expected for such high Mach number shocks because their greater compression ratios enhance the 
acceleration efficiency to the highest energies. 

The complete shock structure is shown in Figure |6|, where we have included along with u x (x) , 
the z -component of flow speed, u z (x) , the angle the local magnetic field makes with the shock 
normal, OBn{%) , and the total magnetic field magnitude, Bt ot (x) . The solid lines show the final 
shock structure and this is compared to the initial structure shown by dotted lines. As noted 
above, both the shape and the overall compression ratio must be modified to obtain a self-consistent 
solution, and this translates to a change from the R-H values in the downstream u z 2 , @Bn2 , and B2 
as indicated by the dotted lines. The sharpness of the subshock in Figure |6| at x = is partially an 
artifact of how we smooth and truncate the flow profile between iterations. As mentioned above, 
we average the predicted profile with the previous one and, since we start with a discontinuous 
shock, some sharpness persists. We also set all predicted values of u x to u x 2 at x > . Despite 
this, the subshock is, in fact, quite sharp as we discuss at the end of this section. 

In Figure ^ we show (solid line) the downstream, shock frame distribution function obtained 
in the smooth shock solution just described. It differs considerably from the test-particle solution 
(dotted line) in that the downstream thermal peak is at a lower energy, the temperature is slightly 
lower, and far fewer thermal particles become accelerated. The straight line is the power-law slope 
expected from the test-particle Fermi solution with r = 5 and matches our test-particle solution 
at energies well above thermal and below where the FEB becomes important. The smooth shock 
solution, however, does not attain the test-particle power-law and remains considerably steeper. 
This can be understood by examining the top panel of Figure [5| where it can be seen that at 
(^feb = — 200 Ao there is enough escaping energy flux to smooth the shock further upstream from 
this point. Particles leave the shock at the FEB before feeling the full compression ratio. 

To complete this example, we show in Figure I distribut ion functions at various x -positions, 
i.e., x = — 50Ao (dotted line), x = — 4Ao (dashed line), x = — 0.5Ao (light solid line), and x = +Ao 
(heavy solid line). At observation points far upstream from the shock, only the unshocked thermal 
peak is present since energetic particles from the shock are not able to diffuse against the background 
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flow to reach the observation point. As the upstream observation point moves toward the shock two 
things happen. First, the highest energy particles begin to show their presence, and second, the 
thermal peak from particles which have not yet crossed the shock begins to shift to lower energy 
(see insert where we have plotted the thermal peaks of the x = — 50Ao and x = — 0.5Ao spectra). 
Since we take A oc r g , the diffusion length increases with energy and higher energy particles from 
the shock are able to stream further upstream than low energy ones, thus as the observation point 
moves toward the shock the spectrum fills in from high energy to low (this property was recognized, 
for test-particle situations, by Baring, Ellison and Jones 1994). The shift in the thermal peak arises 
because spectra calculated in the shock frame shift to lower energy as the bulk flow speed falls as the 
shock is approached. The slowing of the bulk flow in the shock precursor also heats the incoming 
particles somewhat before they encounter the sharp subshock lowering the local Mach number. 
These features are well-defined model predictions that can be tested against observations. 

In Figure |9] we show the upstream scale height for particles of various energies. The ordinate 
is the ratio of the flux at x over the flux at x = for a given energy. As expected, the length scale 
(i.e. distance at which the flux e-folds) is largest for the highest energy particles, and the fluxes 
fall off exponentially with distance from the shock (a property of diffusion against the convecting 
flow). It is also important to note that low energy particles (i.e. the 3 and 10 keV examples) 
can have extremely short upstream precursors. Particle detectors on spacecraft being overtaken by 
interplanetary shocks will see very different time profiles depending on the particle energy sampled 
and may, depending on the time integration of the spectrometer (which is usually long compared 
with typical gyroperiods), see a step function increase in intensity at low energies simultaneously 
with a slow rise in high energy particles. This effect may explain some puzzling aspects of recent 
Ulysses observations which have led to the suggestion that a two-stage acceleration mechanism 
operates for pickup protons (Gloeckler et al. 1994). 

As our final example, we show in Figure |l0| a highly oblique shock ( Oeni = 75° ) with M$i = 
Mai = 10 , d FEB = — 20Ao , rj = 5 , and using large-angle scattering. For these parameters little 
acceleration occurs and the shock profile is nearly discontinuous. Nevertheless, the little smoothing 
evident in the top panel of Figure [l(] is enough to reduce the energy flux from being ~ 20% above 
the far upstream value to a constant value. Because of the inefficient acceleration, few particles, 
carrying very little energy flux, escape at the FEB, so there is no need to adjust the compression 
ratio. The final ratio is r = 3.74 , effectively the R-H value. We have done the same calculation 
with the same parameters except using pitch-angle diffusion and find essentially the same profile 
(which is not shown), although some slight differences do show up in the distribution functions. 



The top two curves in Figure 11 are the distribution function from the LAS example (solid 



line) along with the distribution produced in a shock with the same parameters as that shown 



in Figure 10, only using PAD (dashed line). The main difference between the two cases is that 
the PAD distribution is somewhat smoother than the LAS one. The PAD shock is also somewhat 
less efficient in accelerating particles to the highest energies. In general, however, even at this 
large obliquity, the choice of scattering mode does not play a dominant role in determining the 
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acceleration efficiency. The lower two curves (multiplied by 0.01 for clarity) are test-particle results 
and are similar to the distribution functions produced by the smooth shock, but slightly flatter, 
as expected. No portions of these examples elicit the Fermi, test-particle power-law slope (dotted 
line) but are clearly flattening toward it before the FEB causes the spectra to turn over. 

Finally, we comment on the sharpness of the subshock seen in all of the examples we have 
presented. As mentioned above, part of this is due to the scheme we have for iterating the profile 
and insuring rapid convergence. We average the predicted profile with the previous one, make the 



predicted profile monotonic, and set all predicted values of u x (x) to u x2 for x > . In Figure [12 
we show the same final u x (x) plots as shown in the top left panel of Figure ||, the top panel of 
Figure ||, and the top panel of Figure 10, except here they are all plotted on an expanded linear 



distance scale (solid lines). The dotted lines in Figure 12 are the calculated <v x > obtained from 
distributions produced by the simulation using the solid line shock profiles. In all cases, the mean 
velocities are nearly as sharp as our processed profiles except for statistical fluctuations and the 
truncation we impose at x > . Outside of the range shown, u x and <v x > are indistinguishable 
except for noise. With the possible exception of the low Mach number case, a clear subshock exists 
which is considerably less than Ao in width. While comparing u x and <v x > is somewhat artificial 
since <f x > is calculated using u x , the small broadening of the flow speed around x ~ doesn't 
have an appreciable influence on the momentum and energy fluxes. We have performed simulations 
using <v x > instead of u x (x) and the momentum and energy fluxes remain within ~ 10% of the 
conserved values at all x . For our purposes, since we do not attempt to model scales less than a 
convected thermal ion gyroradius, there is essentially no difference between the u x (x) used in the 
simulation and <v x > . 



4.3. Injection and Acceleration Efficiency 

A quantity that is central to the acceleration problem is the efficiency of the Fermi mechanism. 
It can be defined in a variety of ways: we define the acceleration efficiency, e(> E) , at or behind 
the shock as the downstream energy flux above energy E divided by the incoming energy flux, i.e., 

jr\ _ CP{>E)u x2 + Q esc _ 

where P(>E) is the downstream pressure in particles with energies greater than E . The pressure 
is obtained by taking 2/3 of the energy density in the omni-directional distribution, d J/ dE , and 
we take Q esc to be energy independent, i.e., we assume all of the escaping energy flux is carried 
by the highest energy particles. This last assumption will distort e(> E) somewhat at the highest 
energies since particles of varying energies leave the FEB. 

In Figure [l3| we show e(> E) for our three previous non-linear examples, i.e., curve (a) is for 
the LAS case shown in Figure ||, (b) is for the case shown in Figure [?], and (c) is for the LAS, smooth 
shock case shown in Figure O. As is apparent from the figure, large differences in the efficiency 
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depending on @Bni > the Mach number, r\ , and the distance to the FEB occur. Examples (b) and 
(c) both accelerate particles to above 1000 keV and have similar Mach numbers, but the highly 
oblique shock (c), is much less efficient. Examples (a) and (b) have the same 0Bni = 30° but differ 
considerably in Mach numbers (3 versus 20, respectively), with (b) being more efficient because its 
higher Mach numbers generate a larger compression ratio. At this stage of our work, there appears 
to be no simple way to characterize the injection and acceleration efficiency of oblique shocks, 
but the trends with obliquity and Mach number are quite clear. The flattening at the highest 
energies seen in (a) and (b) is an artifact of assuming that Q esc is carried totally by the highest 
energy particles; if a spread of escape energies were possible, corresponding to a range of d FEB , the 
flattening would be replaced by a quasi-exponential turnover. Note also, that the efficiency curves 



in Figure 13 do not quite converge to unity at low energies due to simulation statistics. 



Another point which should be made concerning Figure [l3| is that we choose to define our 
efficiency as a function of energy. If we do not have an independent source of energetic seed 
particles, all accelerated particles must originate as thermal particles and they will be drawn more 
or less continuously from the thermal population - there will be no clear separation between thermal 
and energetic particles. Our identical treatment of the thermal and non-thermal populations is why 
we have defined our efficiency as a function of energy. However, there are ways to qualitatively 
prescribe efficiencies that describe the overall distribution rather than different particle energies. 



For example, while all three cases in Figure |13| have comparable efficiency at around 0.8 keV, this 
energy does not define the approximate juncture between thermal and non-thermal populations 
for cases (a) and (c), whereas it does for case (b). This juncture is roughly represented by the 
upward kinks at energies 3 keV and 1.7 keV for cases (a) and (c) respectively. Hence, the ratio of 
downstream non-thermal to thermal energy densities for the three cases are roughly (a) 0.2, (b) 0.7 
and (c) 0.2; these numbers could be taken as an alternative measure of acceleration efficiency. 



5. DISCUSSION AND CONCLUSIONS 

A host of observational evidence, both direct and indirect, confirms that collisionless shocks in 
space accelerate particles with high efficiency. Possibly a large fraction of all non-thermal particle 
populations in diffuse regions of space are generated by shocks, making shock acceleration one of 
the most important problems in high energy astrophysics. As a step toward a full understanding of 
shock acceleration, we have developed a model which combines nonlinear particle acceleration and 
diffusion with shock dissipation and non-linear hydrodynamics forming the shock structure. This 
paper presents both the details of our simulation technique and representative acceleration results 
as a prelude to a more comprehensive survey of the parameter space associated with modified, 
oblique shocks. While our model is still incomplete, with simplifying assumptions concerning the 
microphysical processes involved, we believe it is the most realistic current solution of the steady- 
state shock acceleration problem. We include (i) a strongly energy-dependent diffusion coefficient 
[see equation (jj])], which models cross-field diffusion, (ii) the ability to model either large-angle 
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scattering or pitch-angle diffusion, (iii) injection from the thermal background with no additional 
free parameters, (iv) the determination of the self-consistent, average shock structure including 
the dynamic effects of accelerated particles on the thermal shock, (v) the dynamic effects of par- 
ticle escape from finite shocks, and (vi) shock drift and compressional acceleration simultaneously. 
Principal results of this paper include downstream spectra (see Figures ^ and 11), properties of 



upstream populations (see Figures || and ||), and acceleration efficiencies (Figure [131). The Monte 
Carlo simulation does not treat (i) a self-consistent determination of the diffusion coefficient from 
wave-particle interactions, (ii) time-dependent effects, (iii) relativistic particles or flow speeds (iv) 
a cross-shock potential due to charge separation, or (v) geometry other than a plane shock. 

Clearly the most important omission is the self-consistent determination of the diffusion coeffi- 
cient. Our results hinge on the energy-dependent form we have assumed for k , which is motivated 
by previous theoretical analyses (see Jokipii 1987; Giacalone, Burgess and Schwartz 1992; Ellison 
Baring and Jones 1995), and also observational constraints at parallel shocks (see Ellison, Mobius, 
and Paschmann 1990 and Ellison and Reynolds 1991). However, the self-consistent determination 
of k requires knowledge of the microphysics and only plasma simulations where the electric and 
magnetic fields are calculated directly from particle motions (e.g. Quest 1988) can give this informa- 
tion. These simulations are extremely demanding computationally and will not, in the foreseeable 
future, be able to adequately model particle acceleration in shocks to astrophysically important 
energies. The recent work of Jokipii et al. (Jokipii, Giacalone, and Kota 1993; Giacalone, Jokipii, 
and Kota 1994; Jokipii and Jones 1996) has shown that if a coordinate is ignorable (as in one- or 
two-dimensional hybrid simulations), cross- field diffusion effects are suppressed, and since cross- 
field diffusion is an essential part of injection and acceleration in shocks (certainly oblique ones) 
three-dimensional simulations must be used to model shocks. No existing computer is capable of 
running realistic three-dimensional plasma simulations over dynamical ranges of energies appropri- 
ate to astrophysical applications. We emphasize that even though we model plane shocks and our 
scattering operator is a gross simplification of the complex plasma processes taking place in shocks, 
the operator is fully three-dimensional, includes cross-field diffusion, and may well produce more re- 
alistic results than current one- or two-dimension plasma simulations. In fact, comparisons between 
our model and spacecraft observations of highly oblique interplanetary shocks (IPSs) (Baring et al. 
1995) suggest that this is the case. The spacecraft observations clearly show that highly oblique 



shocks inject and accelerate thermal particles, a result we can model accurately (e.g. Figures 11 
and [H|), but one which, to our knowledge, all existing one- and two-dimensional plasma simulations 
fail to show (e.g. Liewer, Goldstein, and Omidi 1993; Liewer, Rath, and Goldstein 1995). 

Virtually all analytic models of nonlinear shock acceleration have been restricted to parallel 
shocks, however, Jones and Kang (1995) have recently extended the cosmic ray diffusion-advection 
equation approach to oblique geometry and produced impressive fits to the Ulysses observations 
mentioned above (e.g. Baring et al. 1995). However, all models based on the diffusion approxi- 
mation (i.e. the requirement that particle speeds be large compared to flow speeds) are limited in 
their ability to treat thermal particles and must use additional free parameters to model injection. 
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For example, the parallel shock model of Berezhko et al. (i.e. Berezhko, Yelshin, and Kseno- 
fontov 1994; Berezhko, Ksenofontov, and Yelshin 1995; Berezhko et al. 1995) uses a source term 
for monoenergetic injection at the gas subshock which is treated as a discontinuity. Here, a small 
fraction e of incoming gas is transferred to cosmic rays, the injected particles instantly obtaining a 
superthermal momentum, pj n j . Both e and pi n j are free parameters and the final results depend 
strongly on them. The main advantage of our model is that the Monte Carlo description is not 
restricted to superthermal particles and injection is treated self-consistently regardless of whether 
or not we assume that all particles obey equation (|l|). Once some such scattering description is 
chosen, both the injection rate and the effective injection momentum are fully determined, as is 
the complete shock structure, by the Monte Carlo solution without any additional parameters such 
as e . In fact, there is no 'injection momentum' in our solution since particles are drawn smoothly 
from the background thermal gas. 

The parameters which determine the injection and acceleration efficiency of shocks are the 
obliquity, 0Bni , the strength of cross-field diffusion, ij , the Mach numbers, M s \ and M A % , and 
the size of the shock system (i.e. |c? FBB | )■ These all influence the shock in complex ways and there 
is no simple relationship between them. In general, we can state that the acceleration efficiency (i.e. 
the fraction of energy flux which ends up in high energy particles) increases with: (i) decreasing 
©Bni j (ii) decreasing r/ (i.e. stronger scattering), (iii) increasing Mach number, and (iv) increasing 
shock size. We have found that the differences in the shock structure and acceleration efficiency 
resulting from using either large-angle scattering (LAS) or pitch-angle diffusion (PAD) are generally 
small (e.g. see Figure |2| for the test-particle regime, and Figures || and 11 for full non-linear results) 



for parameter regime we have investigated, namely v HT <C c, and can be neglected in the Bohm 
diffusion limit. However, the differences between the nonlinear results and the test-particle ones 
are very large (e.g. Figure |7|) except for high obliquities (i.e. Figure 11) or very low Mach numbers 



(Figure B), where the acceleration efficiency is low enough for the thermal gas to dominate the 
non-thermal population dynamically, and the shock profiles are very sharp. 

For the first time, we have been able to calculate the absolute injection and acceleration 
efficiency of non-linear oblique shocks without the use of an ad hoc injection parameter. Our results 
(Figure |l3|) show how large differences in efficiency can occur as parameters change, however, we 
have not yet explored the vast parameter regime oblique geometry opens up. Our next step toward 
a more complete solution of the shock acceleration problem will be a survey intended to quantify the 
differences the various parameters make in the distribution function and overall efficiency. This will 
include determining the effect of varying equation (|l|) and will yield predictions for future spacecraft 
and plasma simulation results. The only way to constrain our r\ parameter is by comparing our 
results with direct observations of shocks or with three-dimensional plasma simulations. Since, 
to our knowledge, no three-dimensional simulation results showing significant acceleration exist, 
we will concentrate on spacecraft observations as they become available. As already mentioned, 
this work has begun with spectral comparisons to Ulysses observations of nearby interplanetary 
shocks (i.e. those not expected to encounter pick-up ions), and our preliminary comparisons with 
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data from Ulysses have already indicated that strong scattering accompanies highly oblique IPSs, 
constraining n to values smaller than about 10 (Baring et al. 1995). Our simulation produces 
particle distributions at different distances upstream and downstream of the subshock, thereby 
providing a wealth of model predictions for testing against observations. 

Once we are successful in constraining our parameters with heliospheric shock observations, 
the next step will be to apply our results to shocks with no in situ particle observations, such as the 
termination shock and supernova remnant blast waves. This is another useful aspect of our model 
and we expect to be able to make predictions for the relative injection and acceleration efficiency 
as a function of ionic composition for both thermal and pickup ions at the termination shock, and 
to calculate how efficiency varies around the rim of a supernova remnant shock. Since relativistic 
particles are produced in these shocks, our model will soon be generalized to include relativistic 
particle energies. 
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A. APPENDIX 



This Appendix describes four technical aspects of the Monte Carlo simulation, namely (1) 
our description of pitch angle diffusion, (2) how particles cross the shock and grid points (actually 
planes of discontinuity of the flow and field profiles) in the simulation, (3) the details of how the 
probability of particle return from beyond the downstream simulation "boundary" is determined, 
and (4) how we derive the form of the conservation of energy flux in equation (^) . 



A.l. Pitch- Angle Diffusion 

To summarize our implementation of pitch-angle diffusion (PAD), which has been given in 
detail in Ellison, Jones, and Reynolds (f990), we simulate small-angle scattering effects by allowing 
the tip of the particle's fluid frame momentum vector p to undergo a random walk on the surface 
of a sphere. If the particle originally had a pitch angle, 9% = arccosjp • B/|p| |B|}, and after a 
time St undergoes a small change in direction of magnitude 59 , its new pitch angle, 9^ , is related 
to the old by 

cos #g = cos 9° cos 59 — sin #° sin 59 cos (f> (Al) 

where cj) is the azimuthal angle of the momentum change 5p measured relative to the plane defined 
by the original momentum p and B . After each scattering, a new phase angle around the magnetic 
field, </>g , is determined from the old phase angle, , by 

sin <b sin 59 



6° + arcsin 



sm9l 



(A2) 



59 is randomly chosen from a uniform distribution between and 59 max , and <j) is randomly 
chosen from a uniform distribution between — it and tt , so that the tip of the momentum vector 
walks randomly over the surface of a sphere of radius p = |p| . 

If the time required to accumulate deflections of the order of 90° is identified with the collision 
time, t c , using a diffusion analysis, the relation between <5# max and the mean free path A was 
shown by Ellison, Jones, and Reynolds (1990) to be 




S9 max = J 6^- (A3) 



where t c = X/v . Pitch angle diffusion is then defined by the regime St <t c . Clearly, using this 
approach implies a magnetic fluctuation correlation length smaller than the particle gyroradius; this 
method then becomes an approximation that is nevertheless still very convenient for implementation 
in Monte Carlo simulations. Note that in the limit of 59 max — > 2ir , this prescription of PAD becomes 
comparable to our scheme for large angle scattering (see Ellison, Jones, and Reynolds 1990). 



Equations ( |A1| ) and ( A2) then can be used to determine the coordinates of the new gyrocenter 

7T 

2 



x — r„ cos 9™ cos f — 5- ) sin Ben 
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y gc = y+ r 6 sin^sin(^-|) (A4) 

IN „ ( jN 7I " 



z gc = z + r g sin 6> B cos ^</> B — J cos ©Bn 

where (x,y,z) is the position of the particle when it scatters. The phase offset of 7r/2 represents the 
differences between position and momentum vector phases. The gyroradius r g remains unchanged 
in the PAD event (true also for large angle scattering) since |B| is fixed at the point of scattering 
and our assumption of elastic scattering leaves the magnitude of the momentum unchanged in the 
fluid frame. However, in contrast to the LAS case where the particle's momentum vector is only 
updated after t c on average, the momentum vector is updated after every 5t for PAD. 



A. 2. Shock or grid zone boundary crossing 

When a particle crosses a grid zone boundary (there is no distinction in our code between 
the shock and any other grid boundary), its orbit is changed because the magnetic field changes 
direction and magnitude. The new values of the particle's pitch angle are obtained from the 
assumption that the momentum in the HT frame remains unchanged at the zone boundary; this 
follows from of the absence of drift electric fields in this frame. This method does not require that 
the magnetic moment be conserved; differences between gyrohelix computations at a flow interface 
and the adiabatic approximation are discussed by Terasawa (1979). 

The detailed calculation (see also Decker 1988; Begelman and Kirk 1990; Ostrowski 1991; 
and Takahara and Terasawa 1991) is as follows (the geometry is illustrated in Figure |i~4| ). The 
component of the old momentum (i.e. the momentum before crossing the grid zone boundary) in 
the y-direction is given by 

Py = p^sin^ (A5) 
and the component along the z'-direction is 

Pz' = P±cos<p° (A6) 

where p° is the component of momentum along Bj, and z 1 is perpendicular to the y-B; plane. Note 
that all momenta in this section are measured in the HT frame. The component of momentum 
along the ^''-direction (i.e. the axis perpendicular to the y-Bj+i plane), is given by 

Pz" = Pz' cos AG Bn - Pb sin A©Bn (A7) 

where A0Bn is the difference in 0Bn across the grid zone boundary. The total momentum perpen- 
dicular to the new magnetic field direction, B; + i, is 



p± = \jp\ +pI 



(A8) 
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and the new momentum parallel to the new magnetic field direction is 

Pb = pgcosAG B n+iVsinAe Bn (A9) 
Finally, the new phase around Bi+i is given by 



£ = arctan-^ (A10) 



A. 3. The Probability of Return Calculation 

The details of how particle return from the far side of the downstream return boundary (DRB) 
is effected are presented here; such return boundaries are found not only in our Monte Carlo 
technique, but also in hybrid plasma simulations (e.g Bennett and Ellison 1995). Assume a uniform 
flow with a component of velocity in the positive x -direction of ii x 2 and assume that particles in 
the local fluid frame are isotropic and of speed v F . Quantities denoted by subscript F are measured 
in this fluid (plasma) frame. The flux of particles crossing a y- z plane that is parallel to, and 
downstream of, the subshock interface is proportional to v xs k > the x -component of particle speed in 
any frame in which the shock is at rest. Here 7J xs k > ( < ) for transmissions to the downstream 
(upstream) side of the DRB. 

The probability that particles return to the DRB after crossing it from the upstream side 
is therefore simply (e.g. see Jones and Ellison 1991) the ratio of the flux of particles moving 
upstream of the DRB to the flux of particles moving to the downstream side of this plane. Clearly 
< w xs k < v F + u x 2 defines downstream crossings of the DRB, while —v F + u x 2 < f xs k < 
prescribes upstream crossings. We confine the discussion to cases with v F > u x 2 , and further 
specialize hereafter to the specific shock rest frame where the flow has zero component of velocity 
in the plane of the DRB. Integrating over the angle of the particle velocity relative to the shock 
normal, or alternatively over v xs k , the probability of return P ret to the DRB for isotropic particles 
of speed v F (in the fluid frame) is 



ret 



/. 





■Vxsk^xsk 

V F +U X 2 



^xsk^xsk 







V F + M X 2 



(All) 



This expression is valid for any shock obliquity and is relativistically correct (Peacock 1981). It 
applies to all v F > u x 2 (for v F < u x 2 , P rc t = ) , and the only requirements for its validity are that 
the particles be isotropic in the local fluid frame and that they do not change speed in the region to 
the right of the return plane. To ensure isotropy, we only apply equation ([All ) after particles have 



scattered at least once in the downstream region when LAS is used, or that particles have diffused 
through 90° in the downstream region when PAD is used. The decision for return, or otherwise, 
is made via a random number generator. 
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Now, consider only those particles which return back across the return plane. Both their pitch 
angle relative to the magnetic field, 6 B = arccos{p-B/|p| |B|} , and their phase around the magnetic 
field, cf) B must be determined in the fluid frame. Determining these requires knowledge of u xs k for 
the returning particles, which can be computed by again noting that the flux of particles returning 
across the return plane (moving in the negative x -direction) is proportional to v xs k • This means 
that the number of particles returning with w xs k between v xs ^ and t> xs k + dv xs ^ is proportional to 
Vxskdvxsk . So, for a particular v F , returning particles are drawn from a distribution such that 

7 —^l / *4k*4k = / dN R = N R (A12) 

(Up - U x2 ) Jv xsk JO 

where N R is a random number uniformly distributed between and 1. Therefore, 

v xs k = VNr{-v f + u x2 ) , (A13) 

and from this the x -component of speed in the plasma frame 

v x f = w xsk - ii X 2 = -Vfv^Vr, - u x2 (1 - \/Nr) (A14) 

can be obtained. Choosing a series of random numbers, N R , between and 1 gives the proper 
distribution of returning particles. 

To determine 9 B and (j) B in the fluid frame, we assume that particles will return distributed 
symmetrically around the x-axis in the fluid frame, a consequence of isotropy. This symmetry 
appears also in our special shock frame (see Figure where the flow is orthogonal to the DRB, 
since phases are preserved in velocity transformations orthogonal to this plane. The velocity vector 
of a returning particle will make an angle X in the fluid frame with the x -axis, given by cos 8 X = 
v xf/v f = (vxsk—u^/vp . If it also has an azimuthal angle about the x-axis of (p x , chosen randomly 
between and 2ir , then from Figure using spherical triangles, the values of the fluid frame pitch 
angle 8 B and phase cp B can be expressed in terms of 9 X , cp x and 0Bn2 • Subsequently, the returning 
particle is placed at the downstream return plane and, using the new fluid frame phase and pitch 
angle (and also a new position for the guiding center), propagated upstream by transforming to the 
de Hoffmann- Teller frame. The statistical prescription in this subsection guarantees that our code 
effectively simulates an infinite region to the downstream side of the probability of return plane. 



A. 4. The energy flux conservation equation 

While the pressure terms in the momentum fluxes in equations (0) and (||) are elementary to 
write down, derivation of the forms for the corresponding terms in the energy flux in equation (|9|) 
involve some subtleties. As outlined Section 3.1, we restrict the analysis to gyrotropic particle 
distributions in local fluid frames in the interest of simplicity; we believe that such an approximation 
is quite good, and certainly better than the assumption of isotropy that is ubiquitous in flux 
conservation equation usage in the literature. 
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In a local fluid frame somewhere in the shock environs, consider coordinate axes oriented so 
that the x-direction is aligned with the local magnetic field, but such that a single rotation about 
the y-axis produces an identical orientation to the system depicted in Figure [j]. In this coordinate 
system, a gyrotropic plasma has a diagonal pressure tensor, namely P xx = Pi and P yy = P zz = P± 
with Pij = otherwise. Pn and P± are components of pressure parallel to and orthogonal to 
the ambient field, and for a thermal plasma are related to analogous temperature components by 
two equations of state. Rotating the axes into alignment with the system in Figure |l| yields a 
non-diagonal pressure tensor V due to mixing of the components: 



V 




COS 0Bn Sin 0Bn 

R | 1 

- Sin 0Bn COS ©Bn 



(A15) 



where 7Z is the rotation matrix. This gives a specific form for the pressure tensor of 
P|| cos 2 G Bn + P_l sin 2 G Bn (P± - Py ) sin e Bn cos 6 Bn 

p - ■ ! o p± o 

(Pl - Pi ) sin e Bn cos Bn Pi cos 2 9 Bn + Pi sin 2 6 Bn 



(A16) 



Since pressure represents the spread of velocities about the mean speed, the pressure tensor is 
invariant under bulk velocity transformations. Hence it follows that equation ( |A16| ) defines the 
pressure in the normal incidence frame, and therefore is directly applicable to the flux conservation 
considerations. The above coordinate rotation therefore yields the relationships 



P|| = P xx + P xz tan 6 B n 
Pi = P xx + P xz cot G Bn 



(A17) 



which define the components of the pressure tensor on the fluid frame in terms of normal incidence 
frame tensor components that can be simply determined using the structure of our simulation. 
At this point it becomes apparent that the assumption of gyrotropy in the fluid frame is indeed 
expedient, since it enables complete specification of the fluid frame pressure tensor using only flux 
quantities measured in the NIF in the x-direction: generalizing from the gyrotropic approximation 
would require construction of coordinate grids in the other directions, thereby complicating the 
simulation immensely, with only marginal gain in physical accuracy. 

The energy flux equation can now be simply constructed from the formalism on p. 56 of Boyd 
and Sanderson (1969). The convective contribution to the energy flux is simply P^u x + P zz it 2 . 
The thermal-type (i.e. velocity spread) term is usually written in the form pk n Tu x /{^ — 1) for 
temperature T , where k B is Boltzmann's constant, and 7 is the ratio of specific heats. For the 
non-thermal application here, prescribing the temperature is inappropriate, so we make use an 
equation of state pk B T = Tr{P} = Pn + 2Pj_ in order to generalize to non-thermal situations by 
converting to pressure formalism. It follows that 



3( 7 - 1) 



iT (P " + 



2P, 



u. 



(7-1) 



-Pxx + 3 P xz (tan B n + 2 cot 9 Bn ) 



(A18) 



-32- 



is the non-convective contribution to the energy flux in the x-direction that results from particle 
pressure. This applies to both the non-relativistic gases considered in this paper (with 7 = 5/3), 
and also more general cases with different equations of state (i.e. 4/3 < 7 < 5/3), including 
relativistic plasmas. Using the magnetic field identity B z /B x = tan@Bn then results in the form 
given in equation (||). 
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Fig. 1. — The geometry of a discontinuous plane shock in the normal incidence frame (NIF), the 
frame where particle distributions and fluxes are output from the simulation. The shock lies in the 
y - z plane and all quantities vary only in x . The geometry of the 'smooth' shock is similar except 
that the magnetic field and flow make many small kinks instead of one large one, i.e. the profiles 
consist of a sequence of connected pieces like this depiction. 

Fig. 2. — Test-particle omni-directional, distribution functions measured several Ao downstream 
from a discontinuous shock in the normal incidence frame. All spectra here and elsewhere are 
normalized to one particle per square cm per second injected far upstream. Each pair of histograms 
has obliquity, Mach numbers and r\ as indicated according to the labels (a), (b) and (c). The solid 
lines are results using large-angle scattering, while the dotted lines are for pitch-angle diffusion. 
At energies where v 3> u x i/cos@B n i applies, all spectra attain the canonical Fermi power-law 
(solid lines of arbitrary normalization). Note that the spectra labeled (c) are multiplied by 0.01 for 
clarity. The similarities for the two modes of scattering are striking. 

Fig. 3. — Flow and flux profiles for quantities in the x -direction for a weak shock with parameters 
Ma.i = 3, M s i = 3, Geni = 30deg, d FEB = — 50Ao , r] = 2, and a self-consistently determined 
compression ratio, r = 2.7. The three left panels show the x -component of the flow speed, the 
flux F^ x (x) of the x-component of momentum, and the energy flux for large-angle scattering. All 
quantities are measured in units of the far upstream values (UpS). The three right panels show the 
same quantities for pitch-angle diffusion. In all cases, four iterations are shown; the first with a light 
solid line, the second with a dashed line, the third with a dot-dashed line, and the fourth with a 
heavy solid line (the flat dotted line indicates upstream values). After four iterations, all quantities 
remain the same except for statistical variations. The momentum and energy fluxes are conserved 
everywhere to within 10%. Note that previous iterations (not shown) were done to determine the 
compression ratio. The different scattering modes produce identical results within statistics. 

Fig. 4. — Distribution functions measured downstream from the shock in the shock (or normal 
incident) frame obtained from the profiles shown in Figure |3|. The turnover at E ~ 100 keV is 
the result of particles leaving at the upstream FEB. Apart from statistics, there are no discernible 
differences between the large-angle scattering and pitch-angle diffusion scattering modes. For com- 
parison, the dotted line shows the test-particle spectrum obtained for the same shock parameters 
and the light solid line is the expected Fermi power-law (arbitrary normalization) for r = 2.7. 

Fig. 5. — The top panel shows the x -component of flow speed versus x . The next two panels show 
the flux i^ x (x) of the x-component of momentum, and the bottom two panels show the energy flux. 
The shock parameters are: M a i = M A i = 20, d PEB = — 200Ao , u\ = 500km sec , n\ = 1cm -3 , 
@Bni = 30deg, and rj = 2 . In all cases, four iterations are shown the first, second, third, and 
fourth iterations are shown by light solid lines, dashed lines, dash-dot lines, and heavy solid lines, 
respectively, and the horizontal scale is logarithmic for x < — 10Ao and linear for x > — Ao ■ After 
four iterations, no further changes occur in the profiles and momentum and energy are conserved 
across the shock once the escaping energy flux (Qesc/-^eni — 0.17) is accounted for. Note that 
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previous iterations were performed to obtain the self-consistent compression ratio, r ~ 5 . 

Fig. 6. — Complete shock structure (solid lines) for the example shown in Figure || again with a 
composite distance scale. The dotted lines show the discontinuous shock structure with the initial 
R-H compression ratio, r = 3.96 . The final shock has a overall compression ratio of r ~ 5 . 

Fig. 7. — The downstream, shock frame distribution functions for the smooth shock obtained in 
Figure [5] (heavy solid line) and a test-particle shock (dotted line) with the same parameters. Note 
the shift of the thermal peak to lower energy which occurs in the smooth shock. The light solid line 
shows the test-particle, power-law slope expected from Fermi acceleration for a shock with r = 5 . 
This is obtained by the discontinuous shock solution before the falloff produced by the FEB, but 
not by the smooth shock solution. 

Fig. 8. — Shock frame distribution functions calculated at x = — 50Ao (dotted line), x = — 4Ao 
(dashed line), x = — 0.5Ao (light solid line), and x = +\q (heavy solid line). The insert shows 
the thermal peaks for the x = — 50Ao and x = — 0.5Ao cases on a linear energy scale. The — 4Ao 
thermal peak lies between these two. 

Fig. 9. — Ratio of the omni-directional flux, dJ/dE, at position x to the flux at x = for a 
given energy versus distance upstream from the shock. These define the scale heights for the several 
energies indicated on the figure. 

Fig. 10. — The three panels show the x -component of the flow speed, the flux F xx (x) of the x- 
component of momentum, and the energy flux for a shock with 0Bni = 75deg , M s \ = M A \ = 10 , 
<^feb = — 20Ao , rj = 5 , and large angle scattering. The profile for the PAD case is essentially the 
same and is not shown. In all panels, the first, second, third, and fourth iterations are shown by 
light solid lines, dashed lines, dash-dot lines, and heavy solid lines, respectively. Comparatively 
little acceleration takes place and the final shock profile is very nearly discontinuous with the R-H 
compression ratio, r = 3.74 . 

Fig. 11. — Distribution functions for a highly oblique shock ( GBni = 75deg ). The top two curves 
are smooth shock solutions using LAS (solid line) and PAD (dashed line). The lower two curves are 
test-particle results. The straight dotted line shows the slope expected from the standard Fermi 
power-law. Unlike our examples at lower obliquities, clear differences exist depending on the type 
of scattering although these are not great enough to produce noticeable differences in the shock 
profile. 

Fig. 12. — Comparison of u x (x)/u\ predicted by the simulation (solid lines) with <v x > (dotted 
lines) plotted with a linear distance scale. The top panel is the example shown in the top left 
panel of Figure ||, the middle panel is the shock shown in Figure ||, and the bottom panel is the 



shock shown in Figure 10. While some sharpness of the shock profile is the result of our smoothing 



procedures, the subshock transition is less than IAq wide in all cases. 
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Fig. 13. — The acceleration efficiency, as defined in equation (19), as a function of particle energy. 
The solid curve (a) is calculated from the LAS spectrum shown in Figure |3] ( GBni = 30° , rj = 2 , 
M A i = M s i =3), the dashed curve (b) is calculated from the smooth shock spectrum shown in 
Figure [7| ( Bn i = 30° , r? = 2, M Al = M sl = 20), and curve (c) is calculated from the LAS, 
smooth shock solution shown in Figure |ll| ( 0Bni = 75° , rj = 5 , M A i = M s i = 10 ). The flattening 
evident at the high energy ends of curves (a) and (b) comes about because to 
be a constant, neglecting the fact that escaping particles leave with a range of energies. 

Fig. 14. — The geometry for calculating the gyroradius and phase of particles crossing a grid zone 
boundary. The boundary occurs at the kink in B . 

Fig. 15. — The velocity vectors of returning particles (filled region in figure, A) are symmetric 
about the x-axis. The vertical dashed line is at zero velocity in the shock frame and the center 
of the v F vectors is displaced by the flow speed, u x 2 . Figure B shows the spherical triangle, abc, 
used to calculate 9 B and (j) B ■ Points a and b lie in the x-z plane, while point c lies off the plane. 
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Fig. 6 
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Fig. 8 
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Fig. 9 




-60 -40 -20 



Distance, x [A ] 



pt_obq_scole.com 



-48- 



Large— Angle Scattering 



Fig. 10 



> 0.5 



©Bn1=75° 

M S1 =M A1 = 10 

T} = 5 



H h 



00 



H 1 1 1 1 1 1 1 1 1 1 h 



FEB 



X 



E 

15 



| 0.5 
o 







H 1 1 1 1 1 1 1 1 1 1 1 1 h 



1 .5 



X 

lT 1 
2 0.5 



-20 -10 

Distance, x [X Q ] 



10 



pt_o b q_d 75_los.com 



-49- 



Fig. 11 
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